{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analyitc solution for CSTR of 2 A -> B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "symbs = t, f, fcA, fcB, IA, IB, k, c1, c2 = sp.symbols('t f phi_A phi_B I_A I_B k c1 c2', real=True, positive=True)\n",
    "symbs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytic(t, f, fcA, fcB, IA, IB, k, c1, c2):\n",
    "    u = sp.sqrt(f*(f + 4*fcA*k))\n",
    "    v = u*sp.tanh(u/2*(t - c1))\n",
    "    return [\n",
    "        (-f + v)/(2*k),\n",
    "        sp.exp(-f*t)*c2 + (f + 2*(fcA + fcB)*k - v)/(2*k)\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exprs = analytic(*symbs)\n",
    "exprs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exprs0 = [expr.subs(t, 0) for expr in exprs]\n",
    "exprs0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = sp.solve([expr - c0 for expr, c0 in zip(exprs0, [IA, IB])], [c1, c2])\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exprs2 = [expr.subs(dict(zip([c1, c2], sol[1]))) for expr in exprs]\n",
    "exprs2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[expr.subs(t, 0).simplify() for expr in exprs2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from chempy import ReactionSystem\n",
    "from chempy.kinetics.ode import get_odesys\n",
    "rsys = ReactionSystem.from_string(\"OH + OH -> H2O2; 'k'\")\n",
    "cstr, extra = get_odesys(rsys, include_params=False, cstr=True)\n",
    "fr, fc = extra['cstr_fr_fc']\n",
    "print(cstr.names, cstr.param_names)\n",
    "cstr.exprs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "args = tout, c0, params = np.linspace(0, .17), {'OH': 2, 'H2O2': 3}, {'k': 5, fc['OH']: 42, fc['H2O2']: 11, fr: 13}\n",
    "res = cstr.integrate(*args)\n",
    "res.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytic_alt(t, f, fcA, fcB, IA, IB, k):\n",
    "    u = sp.sqrt(f*(f + 4*fcA*k))\n",
    "    q = sp.atanh(-u*(sp.sqrt(f) + 2*k*IA)/(f*(f + 4*fcA*k)))\n",
    "    v = u*sp.tanh(u/2*t - q)\n",
    "    w = sp.exp(-f*t)/2/k\n",
    "    y = 2*k*(fcA + fcB)\n",
    "    return [\n",
    "        (-f + v)/(2*k),\n",
    "        w*(sp.exp(f*t)*(f + y - v) - y + 2*k*(IA + IB))\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytic_alt2(t, f, fcA, fcB, IA, IB, k, n):\n",
    "    one_point_five = sp.S(3)/2\n",
    "    a0, b0, x, y = fcA, fcB, IA, IB\n",
    "    Sqrt, Tanh, ArcTanh, E = sp.sqrt, sp.tanh, sp.atanh, sp.E\n",
    "    return [\n",
    "        (-f + Sqrt(f)*Sqrt(f + 8*a0*k)*\n",
    "     Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*t - 2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/2)\n",
    "    )/(4*k),\n",
    "    \n",
    "    (-8*b0*k + 8*b0*E**(f*t)*k + E**(f*t)*f*n - 4*a0*k*n + 4*a0*E**(f*t)*k*n + 4*k*n*x + 8*k*y - \n",
    "    E**(f*t)*Sqrt(f)*Sqrt(f + 8*a0*k)*n*Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*\n",
    "         (t - (2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/(Sqrt(f)*Sqrt(f + 8*a0*k))))\n",
    "        /2))/(8*E**(f*t)*k)\n",
    "    ]\n",
    "#     return [\n",
    "#     (-f + Sqrt(f)*Sqrt(f + 8*a0*k)*\n",
    "#    Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*t - 2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/2)\n",
    "#   )/(4*k),\n",
    "# (E**(f*t)*f - 4*a0*k - 8*b0*k + 4*a0*E**(f*t)*k + 8*b0*E**(f*t)*k + 4*k*x + 8*k*y - \n",
    "#   E**(f*t)*Sqrt(f)*Sqrt(f + 8*a0*k)*Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*\n",
    "#        (t - (2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/(Sqrt(f)*Sqrt(f + 8*a0*k))))\n",
    "#       /2))/(8*E**(f*t)*k)\n",
    "#     ]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = sp.Symbol('n')\n",
    "exprs_alt = analytic_alt2(*symbs[:-2], n)\n",
    "exprs_alt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cses, expr_cse = sp.cse([expr.subs({fcA: sp.Symbol('fr'), fcB: sp.Symbol('fp'), f: sp.Symbol('fv'),\n",
    "                                     IA: sp.Symbol('r'), IB: sp.Symbol('p')}) for expr in exprs_alt])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    '\\n'.join(['%s = %s' % (lhs, rhs) for lhs, rhs in cses] + ['return (\\n    %s\\n)' % str(expr_cse)[1:-1]]).replace(\n",
    "        '3/2', 'three/2').replace(\n",
    "        'exp', 'be.exp').replace(\n",
    "        'sqrt', 'be.sqrt').replace(\n",
    "        'atanh', 'ATANH').replace(\n",
    "        'tanh', 'be.tanh\\n    ').replace(\n",
    "        'ATANH', 'atanh')\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[expr.subs(t, 0).simplify() for expr in exprs_alt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[expr.diff(t).subs(t, 0).simplify() for expr in exprs_alt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(list(rsys.substances))\n",
    "rsys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cstr.names, cstr.param_names)\n",
    "cstr.exprs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "symbs[:-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_cb = sp.lambdify(symbs[:-2] + (n,), exprs_alt)\n",
    "def calc_analytic(xout, y0, p):\n",
    "    return _cb(xout, p[fr], p[fc['OH']], p[fc['H2O2']], y0['OH'], y0['H2O2'], p['k'], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_analytic(result):\n",
    "    ref = calc_analytic(\n",
    "        result.xout,\n",
    "        {k: res.named_dep(k)[0] for k in result.odesys.names},\n",
    "        {k: res.named_param(k) for k in result.odesys.param_names})\n",
    "    return np.array([ref[{'OH': 0, 'H2O2': 1}[k]] for k in result.odesys.names]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yref = get_analytic(res)\n",
    "print(yref.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot()\n",
    "res.plot(y=yref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot(y=res.yout - yref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from chempy.kinetics.integrated import binary_irrev_cstr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_analytic2(result):\n",
    "    ref = binary_irrev_cstr(result.xout, result.named_param('k'), result.named_dep('OH')[0],\n",
    "                           result.named_dep('H2O2')[0], result.named_param(fc['OH']), result.named_param(fc['H2O2']),\n",
    "                           result.named_param(fr))\n",
    "    return np.array([ref[{'OH': 0, 'H2O2': 1}[k]] for k in result.odesys.names]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.plot(y=res.yout - get_analytic2(res))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pseudo reversible reaction (batch reactor)\n",
    "#\n",
    "#     A  +  B   <->   C\n",
    "# t  a0    b0-x     c0 + x\n",
    "import sympy\n",
    "kf, kb, a0, b0, c0, x = sympy.symbols('k_{\\\\rm\\\\ f} k_{\\\\rm\\\\ b} a0 b0 c0 x', real=True, nonnegative=True)\n",
    "rf = kf*a0*(b0-x)\n",
    "rb = kb*(c0+x)\n",
    "dxdt = (rf - rb).expand().factor(x)\n",
    "dxdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integral = sympy.integrate(1/dxdt, (x, 0, x))\n",
    "integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol, = sympy.solve(integral - t, x)\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sympy.ccode(sol.subs({kf: 'kf', kb: 'kb', a0: 'major', b0: 'minor', c0: 'prod'})))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
